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ABSTRACT 

Extending the formalism of Datta, Bharadwaj & Choudhury (2007) for detecting ionized bub- 
bles in redshifted 21 cm maps using a matched-filtering technique, we use different simula- 
tions to analyze the impact of HI fluctuations outside the bubble on the detectability of the 
bubble. In the first three kinds of simulations there is a spherical bubble of comoving radius 
Rfj, the one that we are trying to detect, located at the center, and the neutral hydrogen (HI) 
outside the bubble traces the underlying dark matter distribution. We consider three differ- 
ent possible scenarios of reionization, i.e., (i) there is a single bubble (SB) in the field of 
view (FoV) and the hydrogen neutral fraction is constant outside this bubble (ii) patchy reion- 
ization with many small ionized bubbles in the FoV (PRl) and (iii) many spherical ionized 
bubbles of the same radius Rb (PR2). The centers of the extra bubbles trace the dark matter 
distribution. The fourth kind of simulation uses more realistic maps based on semi-numeric 
modelling (SM) of ionized regions. We make predictions for the currently functioning GMRT 
and a forthcoming instrument, the MWA at a redshift of 6 (corresponding to a observed fre- 
quency 203 MHz) for 1000 hrs observations. We find that for both the SB and PRl scenarios 
the fluctuating IGM restricts bubble detection to size Rb < Q Mpc and Rb < 12 Mpc for 
the GMRT and the MWA respectively, however large be the integration time. These results 
are well explained by analytical predictions. In the PR2 scenario, we find that bubble detec- 
tion is almost impossible for neutral fraction a;Hi < 0.6 because of large uncertainty due to 
the HI fluctuations. Applying the matched-filter technique to the SM scenario, we find that it 
works well even when the targeted ionized bubble is non-spherical due to surrounding bub- 
bles and inhomogeneous recombination. We find that determining the size and positions of 
the bubbles is not limited by the HI fluctuations in the SB and PRl scenario but limited by the 
instrument's angular resolution instead, and this can be done more precisely for larger bubble. 
We also find that for bubble detection the GMRT configuration is somewhat superior to the 
proposed MWA. 
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1 INTRODUCTION 

The epoch of reionization is one of the least known chapters in 
the evolutionary history of the Universe. Quasar absorption spectra 
(Becker et al., 2001; Fan et al., 2002, 2006a) and CMBR observa- 
tions (Spergel et al., 2006; Page et al., 2006; Dunkley et al., 2008) 
together imply that reionization occurred over an extended period 
spanning the redshift range 6 < 2; < 15 (for reviews see Fan, 
Carilli & Keating 2006; Choudhury & Ferrara 2006). It is believed 
that the ionizing radiation from quasars and the stars within galax- 



ies reionize the surrounding neutral intergalactic medium (IGM). 
Ionized bubbles form around these luminous objects, grow and fi- 
nally overlap to completely reionize the Universe. The issue of de- 
tecting these bubbles in radio-interferometric observations of red- 
shifted HI 21 cm radiation has been drawing considerable atten- 
tion. This is motivated by the Giant Metre-Wave Radio Telescope 
(GMRT^; Swamp et al. 1991) which is currently functional and 
several low frequency radio telescopes which are expected to be- 
come functional in the future (eg. MWA^, LOFAR^, 21 CM A 
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PAPERS, VLA extension program'' and SKA^). These are all being 
designed to be sensitive to the HI signal from the epoch of reion- 
ization. 

The detection of individual ionized bubbles would be a direct 
probe of the reionization process. It has been proposed that such ob- 
servations will probe the properties of the ionizing sources and the 
evolution of the surrounding IGM (Wyithe & Loeb, 2004a; Wyithe, 
Loeb & Barnes, 2005; Kohler et al, 2005; Maselli et al., 2007; Al- 
varez & Abel, 2007; Geil & Wyithe, 2008; Wyithe, 2008; Geil et 
al., 2008). Observations of individual ionized bubbles would com- 
plement the study of reionization through the power spectrum of 
HI brightness temperature fluctuations (Zaldarriaga, Furlanetto & 
Hemquist, 2004; Morales & Hewitt, 2004; Bharadwaj & Ali, 2005; 
Ali, Bharadwaj & Pandey, 2005; Sethi, 2005; Datta, Choudhury & 
Bharadwaj, 2007). 

Nearly all the previous work on detecting ionized regions con- 
sider the contrast between the ionized regions and the neutral IGM 
in images of redshifted HI 21 cm radiation. The HI signal is ex- 
pected to be only a small contribution buried deep in the emis- 
sion from other astrophysical sources (foregrounds) and in the sys- 
tem noise. It is a big challenge to detect the signal of an ionized 
bubble from the other contributions that are orders of magnitude 
larger (Shaver et al., 1999; Di Matteo et al., 2002; Oh & Mack, 
2003; Cooray & Furlanetto, 2004; Santos, Cooray & Knox, 2005; 
Gleser et al., 2007; Ali et al., 2008). In an earlier work (Datta, 
Bharadwaj & Choudhury 2007, hereafter referred to as Paper I ) we 
have introduced a matched filter to optimally combine the entire 
signal of an ionized bubble while minimizing the noise and fore- 
ground contributions. This technique uses the visibilities which are 
the fundamental quantity measured in radio-interferometric obser- 
vations. Using visibilities has an advantage over the image based 
techniques because the system noise contribution in different vis- 
ibilities is independent whereas the noise in different pixels of a 
radio-interferometric images is not. 

Paper I presents an analytic framework for predicting the ex- 
pected value and the standard deviation a of the matched filter es- 
timator for the detection of a spherical ionized bubble of comoving 
radius Rb- We identify three different contributions to cr, namely 
foregroimds, system noise and the fluctuations in the HI outside 
the bubble that we are trying to detect. Our analysis shows that 
the matched filter effectively removes the foreground contribution 
so that it falls below the signal. Considering the system noise for 
the GMRT and the MWA we find that a 3 a detection will be pos- 
sible for a bubble of comoving radius Rb > 40 Mpc in 100 hrs 
of observation and Rb > 22 Mpc in 1000 hrs of observation for 
both the instruments. The HI fluctuations, we find, impose a funda- 
mental restriction on bubble detection. Under the assumption that 
the HI outside the ionized bubble traces the dark matter we find 
that it is not possible to detect bubble of size Rb < 8 Mpc and 
Rb < 16 Mpc at the GMRT and MWA respectively. Note that the 
matched filter technique is valid for both, a targeted search around 
QSOs as well as for a blind search in a random direction. 

In this paper we validate the visibility based matched filter 
technique introduced in Paper I through simulations of bubble de- 
tection. Our simulations are capable of handling interferometric ar- 
rays with widely different configurations like the GMRT and the 
MWA , the two instruments that we consider here. As mentioned 



^ http://astro.berkeley.edu/~dbacker/eor/ 
® http://www.cfa.haryard.edu/dawn/ 
http://www.skatelescope.org/ 



earlier, the fluctuations in the HI outside the target bubble impose a 
fundamental restriction for bubble detection. The analytic approach 
of Paper I assumes that the HI outside the bubble traces the dark 
matter. In this paper we carry out simulations that incorporate this 
assumption and use these to assess the impact of HI fluctuations 
for bubble detection. We also use the simulations to determine the 
accuracy to which the GMRT and the MWA will be able to deter- 
mine the size and the position of an ionized bubble, and test if this 
is limited due to the presence of HI fluctuations. In a real situa- 
tion a typical FoV is expected to contain several ionized patches 
besides the one that we are trying to detect. We use simulations to 
assess the impact of HI fluctuations for bubble detection in patchy 
reionization scenarios. 

The outline of the paper is as follows. Section 2 presents a 
brief description of how we simulate 21 -cm maps for three differ- 
ent scenarios of the HI distribution, one where the HI traces the 
dark matter and two with patchy reionization. Subsections 2. 1 and 
2.2 respectively discuss how the simulated maps are converted into 
visibilities and how the matched filter analysis is simulated. We 
present our results in Section 3. Subsections 3.1, 3.2 and 3.3 present 
results for bubble detectability, size determination and position de- 
termination under the assumption that the HI outside the bubble 
traces the dark matter Section 3.4 presents results for bubble de- 
tectability in patchy reionization scenarios. We discuss redshift de- 
pendence of bubble detection in section 4 and present our summary 
in section 5. 

For the GMRT we have used the telescope parameters from 
their website, while for the MWA we use the telescope parame- 
ters from Bowman et al. (2007). The cosmological parameters used 
throughout this paper are Q-m = 0.3, Ht/i^ = 0.022, Us = l.,h = 
0.74, as = 1. 



2 METHOD OF SIMULATION 

We have simulated the detection of the HI signal of an ionized 
bubble whose center is at redshift Zc = 6 which corresponds to 
Vc = 203 MHz. The choice of z value is guided by the fact that we 
expect large ionized regions towards the end of reionization z > 6 
(Wyithe & Loeb, 2004b; Furlanetto et al., 2006). Our aim here is 
to validate the analytic calculations of Paper I and hence the exact 
value of z is not very important. 

We consider four scenarios of reionization for bubble detec- 
tion. In the first three scenarios there is a spherical ionized bubble, 
the one that we are trying to detect, at the center of the FoV. This 
bubble has comoving radius Rb and is embedded in HI that traces 
the dark matter. In the first scenario there is a single bubble in the 
field of view. We refer to this as the SB scenario. In this scenario 
the HI fraction xm is assumed to be uniform outside the bubble. 
The uncertainty due to the HI fluctuations is expected to be lowest 
in this scenario because of the absence of patchiness. This is the 
most optimistic scenario for bubble detection. 

In the next two scenarios, we attempt to quantify the effect of 
patchy reionization (PR) outside the bubble that we are trying to 
detect by introducing many other, possibly overlapping, bubbles in 
the FoV. Unfortunately, there is no obvious way to fix the sizes of 
these bubbles from any theoretical models as they depend crucially 
on the nature of reionization sources and other physical factors. 
In scenario PRl, we assume that the large HII regions which we 
are trying to detect are surrounded by many small ionized regions 
whose sizes are fixed by the following procedure: we assume the 
globally averaged neutral fraction xhi to be ~ 0.5; the reason for 
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Figure 1. This shows HI images on slices through the center of the bubble for the four different scenarios SB, PRl, PR2 and SM (from left to right). In first 
three panels the central, circular dark region of radius Rt = 20 Mpc shows the HII bubble that we are trying to detect. The HI outside this bubble traces the 
dark matter distribution. In the SB scenario(left) the hydrogen neutral fraction is xm = 1 outside the bubble. In the PRl scenario (2nd from left) the extra 
bubbles are all of a fixed comoving radius 6 Mpc. In the PR2 scenario (3rd from left) the extra bubbles have the same comoving radius as the bubble that we 
are trying to detect. In both the PRl and PR2 reionization scenarios the centers of the extra bubbles trace the dark matter distribution and xm = 0.62. In the 
SM scenario (right) the central region up to radius 27 Mpc is fully ionized (marked with solid circle) and beyond that region up to radius 42 Mpc the region 
is partially filled with HI patches (dashed circle). The mean neutral fraction is xjji = 0.5. These simulations are all for the GMRT. 



this choice is that the effects of patchiness would be most promi- 
nent when typically half of the IGM is ionized. Given the value of 
Xm, we ti^ to obtain a reasonable estimate of the size of the back- 
ground bubbles from available models. For example, semi-numeric 
simulations of patchy reionization (Mesinger & Furlanetto, 2007) 
predict that the bubble size distribution peaks around 5 Mpc when 
Xm = 0.61 (see their Fig 6). We thus choose the spherical back- 
ground bubbles to have radii 6 Mpc and compute the number of 
background bubbles by demanding that the resulting neutral frac- 
tion is 0.5. The bubble centres are chosen such that they trace the 
underlying dark matter distribution. At the end, the value of xm 
turns out to be slightly higher 0.62 because of overlap of the bub- 
bles. Note that because of these overlaps, the shapes of the resulting 
ionized regions would not always be perfectly spherical. In this sce- 
nario, we have essentially attempted to capture a situation where 
there are many small, possibly overlapping ionized regions pro- 
duced by galaxies and a few large ionized regions (like the one 
that we are trying to detect) produced by QSOs. 

Since the choice of the background bubble size is not robust 
by any means, we consider a different scenario PR2 where these 
bubbles have the same comoving radius Rb as the bubble that we 
are trying to detect. The centers of these extra bubbles trace the 
dark matter distribution as in PRl. The number of bubbles is fixed 
by the globally averaged xm which we take to be 0.62 same as in 
PRl. The PR2 scenario represents a situation where we predomi- 
nantly have large ionized regions produced either by rare luminous 
sources or through the overlap of several small ionized regions in 
the later stages of reionization. 

A particle-mesh (PM) N-body code was used to simulate the 
dark matter distribution. Our earlier work (Paper I) shows that the 
HI signal of the ionized bubble is largely concentrated at small 
baselines or large angular scales, thus a very high spatial resolu- 
tion is not required. We have used a grid spacing of 2 Mpc for the 
simulations. This is adequate for bubbles in the range 4 < i?;, < 
50 Mpc that we consider The simulations use 256"^ particles on 
a 256^ mesh. For the GMRT a single N-body simulation was cut 
into 8 equal cubes of size 256 Mpc on each side. Considering that 
each cube may be viewed along three different directions, we have 
a total of 24 different realization of the dark matter distribution. 
Each cube corresponds to 18 MHz in frequency and ~ 2° in angle 
which is comparable to the GMRT FoV which has FWHM=1.7° 
at 203MHz. The MWA FoV is much larger (FWHM=13°). Here 



eight independent N-body simulations were used. Viewing these 
along three different directions gives twenty four different realiza- 
tions of the dark matter distribution. Limited computer memory re- 

o 

stricts the simulation size and the angular extent 4 ) is consid- 
erably smaller than the MWA FoV. We do not expect this to affect 
the signal but the contribution from the HI fluctuations outside the 
bubble is possibly underestimated for the MWA. 

The dark matter density contrast S was used to calculate 
the redshifted 21 — cm specific intensity Ii, = Ii,xm{l + 
5) for each grid point of our simulation. Here — 2.5 x 

IO't? (^) (¥) (tIiI)) ™d a^Hi the hydrogen neutral frac- 
tion is inside the ionized bubbles and 1 outside. The simulated 
boxes are transformed to frequency and sky coordinate. Figure 1 
shows the HI image on a slice through the center of the bubble of 
radius Rt = 20 Mpc. The mean neutral fraction xhi is ~ 1 in the 
SB scenario, while it is ~ 0.62 for the two PR simulations shown 
here. 

The three scenarios discussed above consider only spherical 
bubbles, and the only departures from sphericity arise from bub- 
ble overlap. It is important to assess how well our bubble detec- 
tion technique works for non-spherical bubbles, which we do us- 
ing ionization maps produced by the semi-numeric (SM) approach. 
In particular, we use maps obtained by the method of Choudhury, 
Haehnelt & Regan (2008). Essentially, these maps are produced 
by incorporating an excursion-set based technique for identify- 
ing ionized regions given the density distribution and the ionizing 
sources (Zahn et al., 2007; Mesinger & Furlanetto, 2007; Geil & 
Wyithe, 2008). In addition, the method of Choudhury, Haehnelt & 
Regan (2008) incorporate inhomogeneous recombination and self- 
shielding of high-density regions so that it is consistent with the 
"photons-starved" reionization scenario implied by the Lya forest 
data (Bolton & Haehnelt, 2007; Choudhury, Ferrara & Gallerani, 
2008). We use a simulation box of size 270 Mpc with 2000^ par- 
ticles which can resolve collapsed halos as small as ~ IO^Mq. 
The ionization maps are generated at a much lower resolution with 
a grid size of 2.7 Mpc. The box corresponds to 19 MHz in fre- 
quency and ~ 2° in angle comparable to the GMRT FoV. We have 
assigned luminosities to the collapsed halos such that the mean neu- 
tral fraction xhi ~ 0.5. The most massive halo (mass ~ W^^AIq) 
identified in the box is made to coincide with the box centre and 
we assume that it hosts a luminous QSO; its luminosity and age 
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are chosen such that it would produce a spherical HII region of 
comoving size « 27 Mpc in a completely homogeneous neutral 
medium [see, e.g., equation (8) of Geil & Wyithe (2008)]. How- 
ever, the actual ionized region is far from spherical both because of 
the surrounding bubbles from other halos and also because of in- 
homogeneous recombination. We find visually from the maps (see 
the rightmost panel of Figure 1) that the HII region is fully ion- 
ized up to radius w 27 Mpc. Beyond that the region is partially 
filled with neutral patches. This patchy ionized region extends up 
to radius ~ 42 Mpc and then merges with the average IGM. The 
fully ionized region and the region with HI patches are marked with 
two circles. We use this box for GMRT as three independent real- 
izations viewing the box along three different directions. For the 
MWA we need a much larger simulation box which requires sub- 
stantially more computing power, beyond the resources available to 
us at present. Hence we do not consider the MWA for this scenario. 

2.1 Simulating Visibilities 

The quantity measured in radio-interferometric observations is the 
visibility V(U ,v) which is related to the specific intensity pattern 
on the sky lv{0) as 

v(u,v) = j (feA{e)i,.{e)e'"'^" (i) 

Here the baseline U = d/\ denotes the antenna separation d pro- 
jected in the plane perpendicular to the line of sight in units of 
the observing wavelength A, is a two dimensional vector in the 
plane of the sky with origin at the center of the FoV, and A[6) 
is the beam pattern of the individual antenna. For the GMRT this 
can be well approximated by Gaussian A{6) = where 
Qa ~ 0.6 ^FWHM and we use the values 1.7° for Qa at 203 MHz 
corresponding to the redshift z = 6 for the GMRT. The MWA 
beam pattern is expected to be quite complicated, and depends on 
the pointing angle relative to the zenith (Bowman et al., 2007). Our 
analysis largely deals with the beam pattern within 2° of the point- 
ing angle where it is reasonable to approximate the beam as being 
circularly symmetric (Figures 3 and 5 of Bowman et al. 2007 ). We 
approximate the MWA antenna beam pattern as a Gaussian. 

We consider 128 frequency channels across 18MHz band- 
width. The image lv{0) at each channel is multiplied with the 
telescope beam pattern A{Q,v\ The discrete Fourier transform 
(DFT) of the product 1^(9) A{9,i') gives the complex visibili- 
ties V{U, v). The GMRT simulations have baselines in the range 
30.5 < L'^ < 3900 which is adequate to capture the HI signal from 
ionized bubbles which is expected to be confined to small baselines 
U < 1000. 

The visibility recorded in radio-interferometric observations 
is actually a combination of several contributions 

V{U, z/) = S{U, v) + HF{U, z/) + N{U, z/) + F{U, u) . (2) 

The signal 80,1/) from an ionized bubble of comoving radius Rb 
embedded in an uniform HI distribution can be analytically calcu- 
lated (Paper I). The solid curve in Figures 2 and 3 show the ex- 
pected signal for = 20 Mpc. The U extent, frequency extent 
and peak value of the signal scale as R^^, Rb and R^ respectively 
for other values of Rb. Note that S{U,u) is real when the bubble 
is at the center of the FoV. 

The data points shown in Figiu-es 2 and 3 are the real part of 
a few randomly chosen visibiUties determined from the simulation 
oi&Rb = 20 Mpc bubble in the SB scenario. The deviations from 



the analytic predictions are due to the HI fluctuations HF{U, v) ie. 
in the SB scenario the HI outside the bubble traces the dark matter 
fluctuations. Notice that these fluctuations are often so prominent 
that the signal cannot be made out. We expect even larger fluctu- 
ations in the other three scenarios which incorporate patchiness of 
reionization. 

The system noise contribution N{U, v) in each baseline and 

frequency channel is expected to be an independent Gaussian ran- 
dom variable with zero mean ({N) = 0) and variance ^ (iV2)is 
independent of U and Vc. We use (Paper I) 

where C has values 0.53Jy and 54.21Jy for the GMRT and flie 
MWA respectively (Paper I). 

The contribution from astrophysical foregrounds F{U, v) is 
expected to be several order of magnitude stronger than the HI sig- 
nal. The foregrounds are predicted to have a featureless, continuum 
spectra whereas the signal is expected to have a dip at (Figure 
3). This difference holds the promise of allowing us to separate the 
signal from the foregrounds. 

2.2 Simulating Signal Detection 

The signal component S{U, v) in the observed visibilities V{U, v) 
is expected to be buried deep in other contributions many of which 

are orders of magnitude larger. Detecting this is a big challenge. 
For optimal signal detection we consider the estimator (Paper I) 

E = Y,S*f{Ua,Vb)V{Ua,yb) (4) 
a,h 

where Sfitjjv) is a filter which has been constructed to detect 
a particular ionized bubble, V{Ua,Vb) refer to the observed vis- 
ibilities and Ua and Vb refer to the different baselines and fre- 
quency channels in the observation. The filter Sf{U,v) depends on 
{Rf,Zc, Gc] the comoving radius, redshift and angular position of 
the bubble that we are trying to detect. We do not show this explic- 
itly, the values of these parameters will be clear from the context. 

The baselines obtained using DFT in our simulations are uni- 
formly distributed on a plane. In real observations, the baselines 
will have a complicated distribution depending on the antenna 
layout and direction of observation. We incorporate this through 
the normalized baseline distribution function pNiU,^) which is 
defined such that d^Udu pn{U , u) is the fraction of data points 
ie. baselines and frequency channels in the interval d^U dv and 
J d^U J dv pn{U = 1. We use the functional forms of pN 
determined in Paper I for the GMRT and the MWA. 

Using the simulated visibilities, we evaluate the estimator as 

E= {AUf AU J2Sf{Ua,n)V{Ua,Ub)pNiUa,l^b) (5) 

a, 6 

where the sum is now over the baselines and frequency channels in 
the simulation. 

The filter Sf{U, v) (Filter I of Paper I) is defined as 



Sf{U, v) = 




(6) 
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Figure 2. This shows the visibihty signal (real part) from a frequency slice through the center of a spherical ionized bubble of comoving radius 20 Mpc 
embedded in HI. The solid curves show the expected signal assuming that the bubble is embedded in uniformly distributed HI. The data points show the 
visibihties for a few randomly chosen baselines from our simulation of the SB scenario. The difference between the data points and the solid curve is due to 
the fluctuations in the HI outside the bubble. Each panel coiTesponds to a different realization of the simulation. 
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Figure 3. Same as the previous figure except that U is fixed at 110 while the frequency varies , and Ai/ = v — Vc- 



where the first term S{U , i^) is the expected signal of the bubble 
that we are trying to detect. We note that this term is the matched 
filter that gives the maximum signal to noise ratio (SNR). The 
second term involving the Heaviside function Q{x) subtracts out 
any frequency independent component from the frequency range 
Uc — B /2 to Uc + B /2. The latter term is introduced to subtract 
out the foreground contributions. The (v/uc)^ term accounts for 



the fact that Pn{U, i') changes with frequency (equivalently wave- 
length). 

We have used the 24 independent realizations of the simula- 
tion for the first three scenarios to determine the mean (E) and 
the variance {{AE)^) of the estimator. The high computational re- 
quirement restricts us to use just 3 realizations for the SM scenario. 
Only the signal is correlated with the filter, and only this is expected 
to contribute to the mean (E). All the other components are uncor- 
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Figure 4. The estimator E (defined in equation 4) for bubble size Rf, ranging from 4 Mpc to 50 Mpc for tlie GMRT in the SB scenario. It is assumed that 
the filter is exactly matched to the bubble. The left panel shows the analytic predictions for the mean estimator {E) and the 3 — cr error-bars due to the HI 
fluctuations. The solid and the dashed lines in the right panel respectively show the {E) and the 3 — cr envelope determined from the simulations. The data 
points in the right panel show E in the individual realizations. 
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Figure 5. Same as the Figure 4 for the MWA. 



related with the filter and they are expected to contribute only to 
the variance {(Ai?)^). The variance is a sum of three contributions 
(Paper I) 

{{^Ef) = {{^Ef)HF + {{AEf)N + {{AEf)FG ■ (7) 

The simulations give an estimate of {{AE)'^) hf the contribution 
from HI fluctuations. We do not include system noise explicitly in 
our simulations. The noise contribution from a single visibility (eq. 
3) is used to estimate {{AE)'^)n (eq. 19 of Paper I). Under the as- 
sumed foreground model, the foreground contribution {{AE)'^)fg 
is predicted to be smaller than the signal and we do not consider it 
here. 



3 RESULTS 

We first consider the detection of an ionized bubbles and the es- 
timation of its parameters in the SB scenario where there is only 
a single bubble in the FoV. We consider the most optimistic situ- 
ation where the bubble is located in the center. In reality this can 
only be achieved in targeted observations of ionized bubbles around 
luminous QSOs. In a blind search, the bubble in general will be lo- 
cated at some arbitrary position in the FoV, and not the center. It 
has akeady been mentioned that the foregrounds can be removed 



by a suitable choice of the filter. Further, the system noise can, in 
principle, be reduced by increasing the observation time. The HI 
fluctuations outside the bubble impose a fundamental restriction on 
bubble detection. 

3.1 Restriction on bubble detection 

We have carried out simulations for different values of the bubble 
radius Rb chosen uniformly at an interval of 2 Mpc in the range 
4 to 50 Mpc. In each case we consider only the most optimistic 
situation where the bubble radius Rf used in the filter is precisely 
matched to Rb- In reality it is necessary to try filters of different 
radius Rf to determine which gives the best match. 

Figures 4 and 5 shows the results for the GMRT and the MWA 
respectively. We compare the analytic predictions of Paper I (left 
panel) with the prediction of our simulations (right panel). The an- 
alytical predictions for the mean value (E) arising from the sig- 
nal and {(AE)^) due to the HI fluctuations are respectively 
calculated using equations (15) and (22) of Paper I. The signal de- 
pends on the bubble radius Rb and the mean neutral fraction which 
is taken to be xui ~ 1- The uncertainty due to the HI fluctuations 
is calculated using the dark matter power spectrum under the as- 
sumption that the HI traces the dark matter. 
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We find that {E) and \J {{AE)^)^^ determined from the 
simulations is in rough agreement with the analytic predictions. 
The mean {E} is in very good agreement for Rt > 6 Mpc, there is 
a slight discrepancy for smaller bubbles arising from the finite grid 



size (2 Mpc in the simulation). The HI fluctuations y {{AE)'^) 
are somewhat underestimated by the simulations. This is more pro- 
nounced for the MWA where the limited box size of our simula- 
tions results in a FoV which is considerably smaller than the ac- 
tual antennas. We note that the 24 different values of E deter- 
mined fro m the diffe rent realizations of the simulation all lie within 
(E) ± 3-^ {{AEY) determined from the analytic predictions. 

The good agreement between the simulation results and the 
analytical predictions is particularly important because each is 
based on several approximations, many of which differ between the 

two methods. Our results show that the effect of these approxima- 
tions, though present, are well under control. The analytical method 
has the advantage that it is very easy to calculate and can be eval- 
uated very quickly at an extremely low computational cost. Un- 
fortunately, its utility is mainly limited to the SB scenario and it 
cannot be easily applied to an arbitrary PR scenario with a compli- 
cated HI distribution. Simulations, though computationally more 
cumbersome and expensive, are useful in such a situation. It is thus 
important to test that the two methods agree for the SB scenario 
where both of them can be applied. Note that the HI fluctuation 
predicted by the SB scenario sets the lower limit for the HI fluctu- 
ation in any of the PR models. It is expected that patchiness will 
increase the HI fluctuations above the SB predictions. 

It is meaningful to attempt bubble detection at, say 3a" con- 
fidence level, only if {E) > 3^J {{AE)^)hf. The HI fluctua- 
tions overwhelm the signal in a situation where this condition is 
not satisfied, and bubble detection is not possible. In a situation 
where this condition is satisfied, an observed value Eo of the esti- 



mator can be interpreted as a 3cr detection if Eo > 3y {{AE)'^). 
The simulations show that a 3 — a detection is not possible for 
Rb <6 Mpc and Rt < 12 Mpc at the GMRT and MWA respec- 
tively. As noted earlier, the HI fluctuations are somewhat under pre- 
dicted in the simulations and the analytic predictions Ri < 8 Mpc 
and Rb < 16 Mpc respectively, are somewhat larger. 

The limitation on the bubble size Rb that can be detected is 
larger for the MWA as compared to the GMRT. This is because of 
two reasons, the first being the fact that the MWA has a very dense 
sampling of the small baselines where the HI fluctuation are very 
large, and the second being the large FoV. In fact, the baseline dis- 
tribution of the experiment has a significant role in determining the 
quantum of HI fluctuations and thereby determining the lower cut- 
off for bubble detection. Looking for an optimum baseline distribu- 
tion for bubble detection is also an issue which we plan to address 
in future. In a situation where the antenna layout is already in place, 
it may possible to tune the filter to reduce the HI fluctuations. 

We have not considered the effect of peculiar velocities 
(Bharadwaj & All, 2004) in our simulations. From equation (22) of 



Paper I we see that the HI fluctuations scale as y (^[AEY)hf oc 
a/CT, where Ci is the HI multi-frequency angular power spectrum 
(MAPS). The C;s increase by a factor ~ 2 due to peculiar ve- 
locities, whereby \J {{AEY)hf goes up by a factor ~ 1.5. This 
increase does not significantly change our results, and is small com- 
pared to the other uncertainties in the PR models. 

The signal {E) and the HI fiuctuations ^J{{AEy)HF both 



scale as oc xhi, and the lower limit for bubble detection is un- 
changed for smaller neutral fractions. 

3.2 Size Determination 

In this subsection we estimate the accuracy to which it will be pos- 
sible to determine the bubble radius Rb. This, in general, is an un- 
known quantity that has to be determined from the observation by 
trying out filters with different values of i?/. In the matched filter 
technique we expect the predicted SNR (only system noise) ratio 



SNR: 



{E) 



(8) 



NS 



to peak when the filter is exactly matched to the signal ie. Rf ~ 
Rb- The solid line in the right panel of Figure 6 shows this for 
Rb = 10 Mpc. We find that the SNR peaks exactly when the filter 
size Rf = 10 Mpc. We propose that this can be used to obser- 
vationally determine Rb. For varying i?/, we consider the ratio of 

the observed value Eo to the expected system noise ^J {{AE)'^)t^s, 
referring to this as the SNR. The Rf value where this SNR peaks 
gives an estimate of the actual bubble size Rb. The observed SNR 
will differ from the predictions of eq. (8) due to the HI fluctuations 
outside the bubble. These variations will differ from realization to 
realization and this can introduce uncertainties in size estimation. 
We have used the simulations to estimate this. 

The left panel of Figure 6 shows the SNR as a function of Rf 
for 4 different realizations of the simulation for the GMRT with 
bubble size Rb = 10 Mpc. We see that for i?/ < Rb the SNR 
shows a very similar behavior in all the realizations, and it always 
peaks around 10 Mpc as expected. For Rf > Rb the behavior of 
the SNR as a function of Rf shows considerable variation across 
the realizations. In some cases the drop in SNR away from the peak 
is quite rapid whereas in others it is very gradual (for example, 
the dashed-dot-dot curve). In many cases there is an spurious extra 
peak in the SNR at an Rf value that is much larger than Rb. These 
spurious peaks do not pose a problem for size determination as they 
are well separated from Rb and can be easily distinguished from the 
genuine peak. 

The error-bars in the right panel of Figure 6 show the 3 — cr 
fluctuation in the simulated SNR determined from 24 realizations 
of the simulation. Note that the fluctuations at different i?/ are cor- 
related. Although the overall amplitude changes from one reaUza- 
tion to another, the shape of the curve in the vicinity of Rf = Rb 
is nearly invariant across all the realizations. In all of the 24 real- 
izations we can identify a well defined peak at the expected value 
Rf = Rb. 

Figures 7 and 8 show the results for a similar analysis with 
Rb = 20 Mpc for the GMRT and the MWA respectively. It is not 
possible to detect a bubble of size Rb = 10 Mpc with the MWA, 
and hence we do not show this. Here again, we find that for all the 
realizations of the simulations the SNR peaks at Rf = Rb. The 
relative variations in the SNR across the realizations is much less 
for Rb = 20 Mpc as compared to 10 Mpc and there are no spurious 
peaks. Also, for the same bubble size the variations are smaller for 
the GMRT as compared to the MWA. We do not find any spurious 
peaks for Rb = 20 Mpc. 

A point to note is that the mean SNR determined from the 
simulations is somewhat smaller than the analytic predictions, both 
being shown in the right panels of Figures 6, 7 and Figure 8. There 
are a couple of reasons that could account for this namely, (i) the 
bubble in the simulation is not exactly a sphere because of the finite 
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Figure 6. The SNR = {E) / \J ((A_B)2)i.jg for 1000 hrs observation with the GMRT as a function of the filter size Rj for the case where the actual bubble 
size is Rj, = 10 Mpc . The left panel shows 4 different realizations of the simulation. The right panel shows the mean SNR and 3 — cr error-bars calculated 
using 24 reaUzations. The solid line shows the analytical predictions. 




grid size and thus the match between the filter and the signal is 
not perfect even when the sizes are same and (ii) the finite box- 
size imposes a minimum baseline beyond which the signal is not 
represented in the simulation. 

Based on our results we conclude that in the SB scenario for 
the GMRT the accuracy to which the bubble size can be determined 
in our simulations is decided by the resolution 2 Mpc and not by the 
HI fluctuations. In reality the limitation will come from the angular 
resolution of the instrument which sets the limit at 0.5 Mpc for the 
GMRT and 8 Mpc for the MWA. 

The height of the SNR peak depends on the neutral fraction 
and it can be used to observationally determine this. We find that 
the HI fluctuations do not change the position of the peak but intro- 
duce considerable variations in its height even if xm — 1. The HI 
fluctuations restrict the accuracy to which the neutral fraction can 
be estimated, an issue that we propose to address in future work. 



3.3 Determining the position 

In the previous two subsections, we have considered cases where 
the bubble's position is known. Here we assume that the bubble's 
size is known and we estimate the accuracy to which its position 
can be determined in the presence of HI fluctuations. The situation 
considered here is a blind search whereas the former is a targeted 
search centered on a QSO. 

In a real situations it would be necessary to jointly determine 
four parameters the bubble radius Rt,, two angular coordinates 
(&x, dy) and the central frequency Uc from the observation. How- 
ever, to keep the computational requirement under control, in this 
analysis we assume that R^y is known. The bubble is placed at the 
center of the FoV and frequency band, and we estimate how well 
the position can be recovered from the simulation. To determine 
the bubble's position we move the center of the filter to different 
positions and search for a peak in the SNR. To reduce the com- 
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MHz 

Figure 9. The SNR for 1000 hrs GMRT observations for a bubble of size Rt = 10 Mpc located at the center of the filed of view. The filter scans along 
0x , 9y 1 f (top. middle, bottom) to determine the bubble's position. The left panel shows results for 4 reahzations of the SB simulation, the right panels show 
the mean (dashed curve) and 3 — cr error-bars determined from 24 realizations of the simulation and the analytic prediction for the mean (solid curve). 



putational requirement, this is done along one direction at a time, 
keeping the other two directions fixed at the bubble's actual center. 
We have also carried out simulations where the bubble is located 
off-center. We do not explicitly show these results because they are 
exactly the same as when the bubble is at the center except for the 



fact that the value of the peak SNR is lower because of the primary 
beam pattern. 

Figure 9 shows the results for Rt — 10 Mpc for the GMRT. 
The left panel shows results for 4 realizations of the simulation, 
the right panels show the mean and 3 — cr determined from 24 
realizations of the simulations and the analytic prediction for the 




mean value. In all cases a peak is seen at the expected position 
matched with the bubble's actual center. The HI fluctuations pose 
a severe problem for determining the bubble's position as it intro- 
duces considerable fluctuations in the SNR. In some cases these 
fluctuations are comparable to the peak at the bubble's actual posi- 
tion (see the dashed line in the upper left panel). The possibility of 
spurious peaks makes it difficult to reliably determine the bubble's 
position. 



We present the results for Ri, = 20 Mpc in Figures 10 and 
1 1 for the GMRT and MWA respectively. The HI fluctuations do 
not pose a problem for determining the position of such bubbles 
using the GMRT. In all the realizations of the GMRT simulations 
there is a peak at the expected position. The FWHM ~ 40Mpc is 
approximately the same along 6x,0y and Uf and is comparable to 
the separation at which the overlap between the bubble and the filter 
falls to half the maximum value. The HI fluctuations does introduce 
spurious peaks, but these are quite separated from the actual peak 
and have a smaller height. We do not expect these to be of concern 
for position estimation. 



The MWA simulations all show a peak at the expected bubble 
position. The FWHM along (~ 60 Mpc ) is somewhat broader 
than that along u (~ 40 Mpc). The low spatial resolution ~ 8 Mpc 
possibly contributes to increase the FWHM along 6. The HI fluc- 
tuations introduce spurious peaks whose heights are ~ 50 % of the 
height of the actual peak. 



3.4 Bubble Detection in Patchy Reionization 

The SB scenario considered till now is the most optimistic scenario 
in which the HI traces the dark matter. The presence of ionized 
patches other than the one that we are trying to detect is expected 
to increase the contribution from HI fluctuations. We first consider 
the PRl scenario where there are several additional ionized bub- 
bles of radius 6 Mpc in the FoV. Figures 12 & 13 show the mean 
value of the estimator (E) and 3 — ct error-bars as a function of _R/ 
for the GMRT and the MWA respectively. These were estimated 
from 24 different realizations of the simulation, using a filter ex- 
actly matched to the bubble. 

We find that the results are very similar to those for the SB sce- 
nario except that the signal is down by 0.6 due to the lower neutral 
fraction (xm ~ 0.62) in the PR scenarios. Ionized bubbles with 
radius Rt — 8 Mpc and — 12 Mpc or smaller cannot be detected 
by the GMRT and MWA respectively due to the HI fluctuations. 
These limits are similar to those obtained in simulations of the SB 
scenario. 

In the PR2 scenario the FoV contains other ionized bubbles of 
the same size as the bubble that we are trying to detect. We find that 
bubble detection is not possible in such a situation, the HI fluctua- 
tions always overwhelm the signal. This result obviously depends 
on number of other bubbles in the FoV, and this is decided by xm 
which we take to be 0.62. A detection may be possible at higher 
Xm where there would be fewer bubbles in the FoV. 

In the SM scenario, the very large computation time restricts 
us from generating several realizations with central ionized regions 
of different sizes. Hence we are unable to study the restriction on 
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Figure 11. Same as the Figure 9 for Ri, = 20 Mpc for the MWA. 
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bubble detection. We have only three realizations all of which have 
the same ionized region located at the center of the box. Based on 
these we find that the mean estimator (i?) is ~ 30 times larger 
that the standard deviation due to HI fluctuations. The detection 
of a bubble of the size present in our simulation (Figure 1) is not 
restricted by the HI fluctuations. We present size determination re- 
sults in Figure 14. We see that the SNR peaks at i?/ = 42 Mpc 
and not at Rj = 27 Mpc. Recall that in the 21-cm map (Figure 
1) we have visually identified the former as the bubble's outer ra- 
dius which includes several small patchy ionized regions towards 
the periphery and the latter is the inner radius which encloses the 
completely ionized region. We see that the matched filter identifies 
the bubble's outer radius. To study the effect of non-sphericity we 
compare our results in Figure 14 with the predictions for a spherical 
bubble of radius = 42 Mpc embedded in uniform HI with the 
same neutral fraction xhi = 0.5. We find that our results for the 
SM scenario follow the spherical bubble prediction up to a filter 
size Rf = 28 Mpc (marked with a vertical line in Figure 14), cor- 
responding to the bubble's inner radius which encloses a perfectly 
spherical ionized region. Beyond this, and upto the outer radius 
of 42 Mpc, the HI is not fully ionized. There are neutral patches 
which introduce deviations from spherical symmetry and cause the 
SNR to fall below the predictions of a spherical bubble beyond 
Rf = 28 Mpc. The deviations from sphericity also broadens the 
peak in the SNR relative to the predictions for a spherical bubble. 

Our results based on the SM scenario show that the matched 
filter technique works well for bubble detection and for determining 
the bubble's size even when there are deviations from sphericity. 
We obtain good estimates for the extents of both, the completely 



ionized region and the partially ionized region. For the SM sce- 
nario. Figure 15 shows how well the bubble's position can be de- 
termined in a blind search. We have followed the same method as 
described for the SB scenario in subsection 3.3. We see that the 
SNR peaks at the expected position. Further, as the bubble size is 
quite large > 27 Mpc there are no spurious peaks. 



4 REDSHIFT DEPENDENCE 

Results shown so far are all at only one redshift z — 6.1t would be 
interesting and useful to have predictions for higher redshifts. How- 
ever, addressing this issue through direct computations at different 
redshifts would require considerable computation beyond the scope 
of this work. Since we find that the analytic predictions of Paper I 
are in good agreement with the simulations of the SB scenario, we 
use the analytic formalism to predict how different quantities are 
expected to scale with increasing z. 

The redshift dependence of some of the quantities like the sys- 
tem noise, the background 21-cm brightness J,^, and the angular 
and frequency extent of a bubble of fixed comoving radius causes 
the SNR to decrease with increasing z. On the other hand the z de- 
pendence of the neutral fraction, the baseline distribution function 
and the effective antenna collecting area acts to increase the SNR 
at higher redshifts. We find that with increasing z both (E) and 



{(A£')2)hf decrease by nearly the same factor so that the re- 
striction on bubble detection does not change significantly at higher 
redshift in the SB scenario. Assuming that the neutral fraction does 
not change with z, the SNR for bubble detection decreases with 
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increasing redshift, the change depending on the bubble size. For 
example, for the GMRT at 2 = 10 the SNR decreases by a factor 
~ 7 and 6 for bubbles of size Rt — 10 and 20Mpc respectively. 
For the MWA this factor is 3 for both these bubble sizes. We expect 
a similar change in the SNR for the patchy reionization scenarios. 
The drop in SNR is slower for the MWA relative to the GMRT be- 
cause the effective antenna collecting area of the MWA increases 
at higher redshifts. 

The SNR is directly proportional to the global neutral frac- 
tion xui which increases with z. The details of how xm and the 
HI fluctuations change with redshift depends on how reionization 
proceeds with time, an issue beyond the scope of this paper. 



5 SUMMARY 

We have used a visibility-based formalism, introduced in Datta, 
Bharadwaj & Choudhury (2007), to simulate the detection of spher- 
ical HII bubbles in redshifted 21 cm maps through a matched- 
filtering technique. The main aim of this work is to use simulations 
to quantifying the limitations for bubble detection arising from the 
HI fluctuations outside the bubble. We have computed the results 
for two instruments, namely, the GMRT and the upcoming MWA. 
Our main conclusions are as follows: 

• In the case where the HI fluctuations outside the bubble trace 
the dark matter distribution (SB scenario), we find that bubbles with 
radius Rt — 6 Mpc and — 12 Mpc or smaller cannot be detected 
by the GMRT and MWA respectively due to the HI fluctuations. 
Note that this limitation is fundamental to the observations and can- 
not be improved upon by increasing integration time. 

• For targeted observations of ionized bubbles, the bubble size 
can be determined to an accuracy limited by the instrument's reso- 
lution; we find that HI fluctuations do not play any significant role. 
However, the HI fluctuations can restrict the accuracy to which the 
neutral fraction can be estimated. In addition, we find that deter- 
mining the position of the bubble in a blind search could be quite 
difficult for small (~ 10 Mpc) bubbles as the HI fluctuations intro- 
duce large fluctuations in the signal; for larger bubbles the accuracy 
is determined by the instrument's resolution. 

• In a scenario of patchy reionization where the targeted HII re- 
gion is surrounded by many small ionized regions of size ~ 6 Mpc 
(PRl scenario), the lower limit for bubble detection is similar to 
that in the SB scenario. Thus the assumption that the HI traces the 
dark matter gives a reasonable estimate of the contribution from HI 
fluctuations if the background ionized bubbles are small ~ 6 Mpc. 
However, the situation is quite different when the surrounding bub- 
bles as of similar size as the targeted bubble (PR2 scenario). The 
large HI fluctuations do not permit bubble detection for a neutral 
fraction xm < 0.6. Thus for xm = 0.6 or lower, bubble detection 
is possible only if the other ionized regions in the FoV are much 
smaller than the bubble that we are trying to detect. 

• The matched filter technique works well for more realistic 
cases based on the semi-numeric modelling of ionized regions 
(Choudhury, Haehnelt & Regan, 2008). Here the bubbles are sub- 
stantially non-spherical because of surrounding bubbles and inho- 
mogeneous recombination. Our method gives a good estimate of 
the size of both the fully ionized and the partially ionized regions 
in the bubble. 

To put our conclusions in an overall perspective, let us con- 
sider an ionized bubble around a luminous QSO at 2; > 6. We ex- 
pect Rb > 30Mpc from studies of QSO absorption spectra (Wyithe 
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Figure 12. The mean (E) and 3 — cr eiTor-bars of the estimator as a function 
of Rf for the GMRT estimated from 24 different realizations of the PRl 
scenario. In all cases the filter is exactly matched to the bubble. 
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Figure 13. Same as the Figure 12 for the MWA. 

& Loeb 2004b; Mesinger & Haiman 2004). It has also been pointed 
out that these bubbles may survive as large "gray fossils" a long 
time after the source has shut down (Furlanetto et al. 2008). It will 
be possible to detect such bubbles only if the background bubbles 
are smaller, say, < 30 Mpc. We find from models of Mesinger & 
Furlanetto (2007) that the typical sizes of ionized regions when 
xhi ~ 0.3(0.1) is ~ 20(70)Mpc. Though these values could be 
highly model-dependent, it still gives us an idea that the bubbles 
around the luminous QSO would be detectable even in a highly 
ionized IGM with, say, xhi ~ 0.3. If the size of the targeted bub- 
ble is larger, then this constraint is less severe. This gives a realistic 
hope of detecting these bubbles at z > 6 with near-future facilities. 

A caveat underlying a large part of our analysis is the assump- 
tion that the bubbles under consideration are perfectly spherical. 
This is note the case in reality. For example, non-isotropic emis- 
sion from the sources (QSOs), density fluctuations in the IGM and 
radiative transfer effects would distort the shape of the bubble. The 
semi-numeric simulations (SM scenario) incorporate some of these 
effects and give an estimate of the impact of the deviations from 
sphericity on bubble detection. This is an important issue which 
we plan to address in more detail in future. In addition, the finite 
light travel time gives rise to an apparent non-sphericity even if the 
physical shape is spherical (Wyithe & Loeb, 2004a; Yu, 2005). This 
effect can, in principle, be estimated analytically and incorporated 
in the filter. We plan to address this effect in a separate publication. 
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Figure 15. Same as Figure 9 for the SM scenario for the GMRT. The x-axis shows the comoving distance of the filter position from the center of the box. The 
three curves respectively show results for a search along three Ox, 9y and u axes. 
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Figure 14. Same as Figure 6 for the SM scenario for the GMRT. The dot- 
ted, dashed dotted and dashed lines show results for thi'ee different real- 
izations. To show the effect of non-sphericity, we compare these results 
with predictions for a spherical bubble of sized Ri, == 42 Mpc embed- 
ded in uniform HI with neutral fraction xm = 0.5(solid line). The vertical 
line al Rf = 28 Mpc shows the radius up to which the bubble is fully 
ionized and the SNR follows the spherical predictions. The SNR peaks at 
Rf = 42 Mpc marked by another vertical line. 
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